Chapter 4 Community composition

load("data/data.Rdata")

4.1 Taxonomy overview

4.1.1 Stacked barplot

# Merge data frames based on sample
merged_data<-genome_counts_filt %>%
  mutate_at(vars(-genome),~./sum(.)) %>% #apply TSS normalisation
  pivot_longer(-genome, names_to = "sample", values_to = "count") %>% #reduce to minimum number of columns
  left_join(., genome_metadata, by = join_by(genome == genome)) %>% #append genome metadata
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>% #append sample metadata
  filter(count > 0) #filter 0 counts

# Create an interaction variable for time_point and sample
merged_data$interaction_var <- interaction(merged_data$sample, merged_data$time_point)

ggplot(merged_data, aes(x=interaction_var,y=count, fill=phylum, group=phylum)) + #grouping enables keeping the same sorting of taxonomic units
    geom_bar(stat="identity", colour="white", linewidth=0.1) + #plot stacked bars with white borders
    scale_fill_manual(values=phylum_colors) +
    facet_nested(. ~ type,  scales="free") + #facet per day and treatment
    guides(fill = guide_legend(ncol = 1)) +
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),
          axis.title.x = element_blank(),
          panel.background = element_blank(),
          panel.border = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          axis.line = element_line(linewidth = 0.5, linetype = "solid", colour = "black")) +
   labs(fill="Phylum",y = "Relative abundance",x="Samples")+
  scale_x_discrete(labels = function(x) gsub(".*\\.", "", x)) +
  facet_wrap(~type, scales = "free", labeller = as_labeller(function(label) gsub(".*\\.", "", label))) #only show time_point label

4.1.2 Phylum relative abundances

phylum_summary <- genome_counts_filt %>%
  mutate_at(vars(-genome),~./sum(.)) %>% #apply TSS normalisation
  pivot_longer(-genome, names_to = "sample", values_to = "count") %>%
  left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
  left_join(genome_metadata, by = join_by(genome == genome)) %>%
  group_by(sample,phylum) %>%
  summarise(relabun=sum(count))

phylum_summary %>%
    group_by(phylum) %>%
    summarise(mean=mean(relabun, na.rm=TRUE),sd=sd(relabun, na.rm=TRUE)) %>%
    arrange(-mean) %>%
    tt()

phylum	mean	sd
p__Bacteroidota	0.380806816	0.202118047
p__Bacillota_A	0.253813952	0.163904042
p__Bacillota	0.118079729	0.146679823
p__Pseudomonadota	0.096901291	0.160581847
p__Campylobacterota	0.053050348	0.092998130
p__Verrucomicrobiota	0.029867960	0.073128363
p__Desulfobacterota	0.023580475	0.036482484
p__Chlamydiota	0.010572704	0.059690269
p__Fusobacteriota	0.010482132	0.028311738
p__Cyanobacteriota	0.009206465	0.016492133
p__Bacillota_C	0.004713164	0.006645703
p__Spirochaetota	0.004009680	0.012308028
p__Bacillota_B	0.002465465	0.004927667
p__Actinomycetota	0.001235741	0.006346852
p__Elusimicrobiota	0.001214079	0.006501752

phylum_arrange <- phylum_summary %>%
    group_by(phylum) %>%
    summarise(mean=mean(relabun)) %>%
    arrange(-mean) %>%
    select(phylum) %>%
    pull()

phylum_summary %>%
    filter(phylum %in% phylum_arrange) %>%
    mutate(phylum=factor(phylum,levels=rev(phylum_arrange))) %>%
    ggplot(aes(x=relabun, y=phylum, group=phylum, color=phylum)) +
        scale_color_manual(values=phylum_colors[rev(phylum_arrange)]) +
        geom_jitter(alpha=0.5) + 
        theme_minimal() + 
        theme(legend.position="none") +
        labs(y="Phylum",x="Relative abundance")

4.2 Taxonomy boxplot

4.2.1 Family

family_summary <- genome_counts_filt %>%
  mutate_at(vars(-genome),~./sum(.)) %>% #apply TSS normalisation
  pivot_longer(-genome, names_to = "sample", values_to = "count") %>% #reduce to minimum number of columns
  left_join(sample_metadata, by = join_by(sample == Tube_code)) %>% #append sample metadata
  left_join(., genome_metadata, by = join_by(genome == genome)) %>% #append genome metadata
  group_by(sample,family) %>%
  summarise(relabun=sum(count))

family_summary %>%
    group_by(family) %>%
    summarise(mean=mean(relabun, na.rm=TRUE),sd=sd(relabun, na.rm=TRUE)) %>%
    arrange(-mean) %>%
    tt()

family	mean	sd
f__Bacteroidaceae	2.215146e-01	0.1392048812
f__Lachnospiraceae	1.406192e-01	0.1085410432
f__Tannerellaceae	1.028077e-01	0.0799387840
f__Helicobacteraceae	5.258990e-02	0.0925721039
f__Mycoplasmoidaceae	3.694491e-02	0.0756423660
f__Erysipelotrichaceae	3.502391e-02	0.0451216980
f__UBA3700	3.418595e-02	0.0557485850
f__Marinifilaceae	2.774048e-02	0.0271999075
f__	2.772765e-02	0.0838376915
f__DTU072	2.700317e-02	0.0529029477
f__Rikenellaceae	2.696900e-02	0.0464477971
f__Enterobacteriaceae	2.691677e-02	0.0918386902
f__Coprobacillaceae	2.551773e-02	0.0892179821
f__Desulfovibrionaceae	2.358047e-02	0.0364824837
f__Ruminococcaceae	1.859938e-02	0.0429483670
f__LL51	1.759437e-02	0.0687507595
f__Rhizobiaceae	1.596807e-02	0.0770633474
f__UBA3830	1.578817e-02	0.0454016753
f__Akkermansiaceae	1.227359e-02	0.0316838689
f__Chlamydiaceae	1.057270e-02	0.0596902690
f__Fusobacteriaceae	1.048213e-02	0.0283117384
f__CAG-239	9.141331e-03	0.0150992357
f__Enterococcaceae	8.420523e-03	0.0466580889
f__Gastranaerophilaceae	7.728770e-03	0.0144404459
f__Oscillospiraceae	6.643218e-03	0.0078105983
f__UBA1997	6.378408e-03	0.0309668631
f__Streptococcaceae	6.366441e-03	0.0342174908
f__UBA1242	4.673359e-03	0.0153730061
f__Brevinemataceae	4.009680e-03	0.0123080282
f__Acutalibacteraceae	3.374550e-03	0.0109560112
f__UBA660	3.196511e-03	0.0139558713
f__Clostridiaceae	2.842205e-03	0.0171339212
f__RUG11792	2.817729e-03	0.0251127757
f__Peptococcaceae	2.465465e-03	0.0049276669
f__MGBC116941	2.160169e-03	0.0093827415
f__Acidaminococcaceae	1.943804e-03	0.0050321931
f__CAG-508	1.807978e-03	0.0064175379
f__Anaerovoracaceae	1.569531e-03	0.0036471440
f__Moraxellaceae	1.485415e-03	0.0097446552
f__RUG14156	1.477695e-03	0.0045847831
f__Staphylococcaceae	1.361709e-03	0.0050863070
f__Elusimicrobiaceae	1.214079e-03	0.0065017516
f__CAG-288	9.490865e-04	0.0060198775
f__Anaerotignaceae	8.989745e-04	0.0040469504
f__CALVMC01	7.516846e-04	0.0043609744
f__Eggerthellaceae	6.407882e-04	0.0021266467
f__Massilibacillaceae	6.235073e-04	0.0016350480
f__Mycobacteriaceae	5.949530e-04	0.0060400054
f__UBA1820	4.736030e-04	0.0013057353
f__Arcobacteraceae	4.604457e-04	0.0050324221
f__CAG-274	4.519746e-04	0.0022028477
f__Burkholderiaceae_C	3.699430e-04	0.0048092594
f__Muribaculaceae	3.617894e-04	0.0009776409
f__UBA932	3.419549e-04	0.0011583178
f__Hepatoplasmataceae	2.989107e-04	0.0038858387
f__Rhodobacteraceae	2.959092e-04	0.0038468195
f__Weeksellaceae	2.771627e-04	0.0031476408
f__Eubacteriaceae	1.646822e-04	0.0006729067
f__Sphingobacteriaceae	1.505775e-04	0.0012460017
f__Devosiaceae	1.489995e-04	0.0015094318
f__Pumilibacteraceae	1.277418e-04	0.0007646754
f__WRAU01	9.603359e-05	0.0012484367
f__Peptostreptococcaceae	2.287338e-05	0.0002973539

family_arrange <- family_summary %>%
    group_by(family) %>%
    summarise(mean=sum(relabun)) %>%
    arrange(-mean) %>%
    select(family) %>%
    pull()

family_summary %>%
    left_join(genome_metadata %>% select(family,phylum) %>% unique(),by=join_by(family==family)) %>%
    left_join(sample_metadata,by=join_by(sample==Tube_code)) %>%
    filter(family %in% family_arrange[1:20]) %>%
    mutate(family=factor(family,levels=rev(family_arrange[1:20]))) %>%
    filter(relabun > 0) %>%
    ggplot(aes(x=relabun, y=family, group=family, color=phylum)) +
        scale_color_manual(values=phylum_colors[-8]) +
        geom_jitter(alpha=0.5) + 
        facet_grid(.~type)+
        theme_minimal() + 
        labs(y="Family", x="Relative abundance", color="Phylum")

4.2.2 Genus

genus_summary <- genome_counts_filt %>%
  mutate_at(vars(-genome),~./sum(.)) %>% #apply TSS nornalisation
  pivot_longer(-genome, names_to = "sample", values_to = "count") %>% #reduce to minimum number of columns
  left_join(sample_metadata, by = join_by(sample == Tube_code)) %>% #append sample metadata
  left_join(genome_metadata, by = join_by(genome == genome)) %>% #append genome metadata
  group_by(sample,genus) %>%
  summarise(relabun=sum(count)) %>%
  filter(genus != "g__")

genus_summary %>%
    group_by(genus) %>%
    summarise(mean=mean(relabun, na.rm=TRUE),sd=sd(relabun, na.rm=TRUE)) %>%
    arrange(-mean) %>%
    tt()

genus	mean	sd
g__Bacteroides	1.352836e-01	0.0923254347
g__Parabacteroides	9.681305e-02	0.0801790006
g__Phocaeicola	6.935126e-02	0.0794715561
g__Mycoplasmoides	3.051289e-02	0.0753256189
g__Helicobacter_J	3.009004e-02	0.0594808174
g__Odoribacter	2.552307e-02	0.0267570773
g__Roseburia	2.384718e-02	0.0559439989
g__NHYM01	2.249986e-02	0.0796951054
g__Alistipes	2.208351e-02	0.0283591693
g__Coprobacillus	2.013859e-02	0.0878851683
g__Agrobacterium	1.596807e-02	0.0770633474
g__Akkermansia	1.227359e-02	0.0316838689
g__Fusobacterium_A	1.038903e-02	0.0283170720
g__Kineothrix	8.801109e-03	0.0409037874
g__Proteus	8.657984e-03	0.0681831315
g__Dielma	8.578108e-03	0.0090836958
g__CAG-95	8.109350e-03	0.0205019603
g__JAAYNV01	7.994602e-03	0.0195404994
g__UBA866	7.260325e-03	0.0293499698
g__Desulfovibrio	7.018155e-03	0.0211479806
g__Enterococcus	7.001300e-03	0.0456600395
g__Ureaplasma	6.432013e-03	0.0138012183
g__Lactococcus	6.366441e-03	0.0342174908
g__Lacrimispora	6.064761e-03	0.0098190433
g__Parabacteroides_B	5.994680e-03	0.0100174097
g__CALXRO01	5.765728e-03	0.0308524398
g__Citrobacter	5.717035e-03	0.0334547978
g__NSJ-61	5.560111e-03	0.0198839325
g__Breznakia	5.504528e-03	0.0237016548
g__Clostridium_AQ	5.326190e-03	0.0121695008
g__Salmonella	5.133583e-03	0.0184556360
g__Bilophila	4.983840e-03	0.0088456017
g__Hungatella_A	4.811802e-03	0.0095549574
g__MGBC136627	4.364796e-03	0.0163003355
g__Escherichia	4.188365e-03	0.0266100578
g__UMGS1251	4.159842e-03	0.0072716746
g__Hungatella	4.116388e-03	0.0190788279
g__Clostridium_Q	4.011996e-03	0.0052127439
g__Brevinema	4.009680e-03	0.0123080282
g__Thomasclavelia	3.902580e-03	0.0109042017
g__Scatousia	3.649941e-03	0.0102727788
g__Enterocloster	3.620384e-03	0.0047339282
g__Mailhella	3.612079e-03	0.0102470769
g__Copromonas	3.599368e-03	0.0050180473
g__Ventrimonas	3.513355e-03	0.0071233140
g__Caccovivens	3.343316e-03	0.0122194941
g__Lawsonia	3.289015e-03	0.0117180962
g__Fournierella	3.217226e-03	0.0062309921
g__Limenecus	3.163279e-03	0.0065731820
g__MGBC133411	3.042418e-03	0.0074576633
g__Mucinivorans	2.900095e-03	0.0373193950
g__Sarcina	2.842205e-03	0.0171339212
g__Acetatifactor	2.738138e-03	0.0058449483
g__Eisenbergiella	2.697412e-03	0.0068331645
g__Bacteroides_G	2.682722e-03	0.0346150378
g__CAJLXD01	2.633818e-03	0.0087495783
g__Blautia	2.558708e-03	0.0061407647
g__C-19	2.275515e-03	0.0048691968
g__Butyricimonas	2.217402e-03	0.0050081070
g__Velocimicrobium	2.201224e-03	0.0066764022
g__CAZU01	2.196385e-03	0.0065944622
g__MGBC116941	2.160169e-03	0.0093827415
g__Intestinimonas	2.081898e-03	0.0039382613
g__Negativibacillus	2.069077e-03	0.0055137501
g__Rikenella	1.985389e-03	0.0037067264
g__Phascolarctobacterium	1.943804e-03	0.0050321931
g__RGIG6463	1.789843e-03	0.0039682804
g__JALFVM01	1.742360e-03	0.0038623852
g__Oscillibacter	1.510459e-03	0.0024992661
g__Pseudoflavonifractor	1.502976e-03	0.0027706264
g__Acinetobacter	1.485415e-03	0.0097446552
g__Citrobacter_A	1.392923e-03	0.0060348874
g__Staphylococcus	1.361709e-03	0.0050863070
g__RGIG4733	1.303790e-03	0.0040480788
g__UBA1436	1.214079e-03	0.0065017516
g__Lachnotalea	1.208197e-03	0.0049169585
g__Ruthenibacterium	1.202021e-03	0.0053632179
g__14-2	1.184643e-03	0.0096299356
g__Beduini	1.173950e-03	0.0025142208
g__Scatocola	1.121193e-03	0.0044975564
g__Faecisoma	1.085585e-03	0.0055596648
g__Enterococcus_A	1.083991e-03	0.0099150522
g__Faecimonas	9.880526e-04	0.0083079244
g__RGIG9287	9.638826e-04	0.0093228400
g__CAG-345	9.490865e-04	0.0060198775
g__Blautia_A	9.208028e-04	0.0029082304
g__RGIG1896	8.343610e-04	0.0062772562
g__CAG-269	7.982448e-04	0.0047176841
g__Marseille-P3106	7.938268e-04	0.0017618833
g__WRHT01	6.429563e-04	0.0026979819
g__Eggerthella	6.407882e-04	0.0021266467
g__IOR16	6.398942e-04	0.0020651222
g__Anaerotruncus	6.265445e-04	0.0016948517
g__RUG14156	6.151780e-04	0.0022147136
g__CHH4-2	6.145042e-04	0.0019997615
g__Corynebacterium	5.949530e-04	0.0060400054
g__Serratia_A	5.860616e-04	0.0076188009
g__CAG-56	4.915432e-04	0.0016341973
g__Merdimorpha	4.736030e-04	0.0013057353
g__MGBC140009	4.679334e-04	0.0024067320
g__CALURL01	4.635287e-04	0.0016737486
g__Aliarcobacter	4.604457e-04	0.0050324221
g__Scatenecus	4.520215e-04	0.0019760876
g__RGIG8482	4.399064e-04	0.0029753981
g__Enterobacter	4.073437e-04	0.0041317730
g__Klebsiella	4.054439e-04	0.0048910857
g__Caccenecus	3.941198e-04	0.0017802371
g__Alcaligenes	3.699430e-04	0.0048092594
g__Plesiomonas	3.633249e-04	0.0027105056
g__HGM05232	3.617894e-04	0.0009776409
g__JAHHSE01	3.592120e-04	0.0014900895
g__Egerieousia	3.419549e-04	0.0011583178
g__Emergencia	3.413630e-04	0.0017450341
g__Enterococcus_B	3.352316e-04	0.0022266910
g__Stoquefichus	3.026072e-04	0.0020503969
g__Hepatoplasma	2.989107e-04	0.0038858387
g__Paracoccus	2.959092e-04	0.0038468195
g__Moheibacter	2.771627e-04	0.0031476408
g__Scatomorpha	2.641015e-04	0.0010184339
g__UBA7185	2.434328e-04	0.0014558191
g__Eubacterium	1.646822e-04	0.0006729067
g__Sphingobacterium	1.505775e-04	0.0012460017
g__Devosia	1.489995e-04	0.0015094318
g__Anaerosporobacter	1.454112e-04	0.0012747262
g__Caccomorpha	1.383122e-04	0.0010540604
g__UBA2658	1.307451e-04	0.0007204965
g__Protoclostridium	1.277418e-04	0.0007646754
g__Angelakisella	1.268479e-04	0.0009221501
g__Cetobacterium_A	9.310217e-05	0.0008718575
g__Rahnella	6.470705e-05	0.0008411917
g__Peptostreptococcus	2.287338e-05	0.0002973539

genus_arrange <- genus_summary %>%
    group_by(genus) %>%
    summarise(mean=sum(relabun)) %>%
    filter(genus != "g__")%>%
    arrange(-mean) %>%
    select(genus) %>%
    mutate(genus= sub("^g__", "", genus)) %>%
    pull()

genus_summary %>%
    left_join(genome_metadata %>% select(genus,phylum) %>% unique(),by=join_by(genus==genus)) %>%
    left_join(sample_metadata,by=join_by(sample==Tube_code)) %>%
    mutate(genus= sub("^g__", "", genus)) %>%
    filter(genus %in% genus_arrange[1:20]) %>%
    mutate(genus=factor(genus,levels=rev(genus_arrange[1:20]))) %>%
    filter(relabun > 0) %>%
    ggplot(aes(x=relabun, y=genus, group=genus, color=phylum)) +
        scale_color_manual(values=phylum_colors[-c(3,4,6,8)]) +
        geom_jitter(alpha=0.5) + 
        facet_grid(.~type)+
        theme_minimal() + 
        labs(y="Family", x="Relative abundance", color="Phylum")

4.3 Alpha diversity

# Calculate Hill numbers
richness <- genome_counts_filt %>%
  column_to_rownames(var = "genome") %>%
  dplyr::select(where(~ !all(. == 0))) %>%
  hilldiv(., q = 0) %>%
  t() %>%
  as.data.frame() %>%
  dplyr::rename(richness = 1) %>%
  rownames_to_column(var = "sample")

neutral <- genome_counts_filt %>%
  column_to_rownames(var = "genome") %>%
  dplyr::select(where(~ !all(. == 0))) %>%
  hilldiv(., q = 1) %>%
  t() %>%
  as.data.frame() %>%
  dplyr::rename(neutral = 1) %>%
  rownames_to_column(var = "sample")

phylogenetic <- genome_counts_filt %>%
  column_to_rownames(var = "genome") %>%
  dplyr::select(where(~ !all(. == 0))) %>%
  hilldiv(., q = 1, tree = genome_tree) %>%
  t() %>%
  as.data.frame() %>%
  dplyr::rename(phylogenetic = 1) %>%
  rownames_to_column(var = "sample")

# Aggregate basal GIFT into elements
dist <- genome_gifts %>%
  to.elements(., GIFT_db) %>%
  traits2dist(., method = "gower")

functional <- genome_counts_filt %>%
  column_to_rownames(var = "genome") %>%
  dplyr::select(where(~ !all(. == 0))) %>%
  hilldiv(., q = 1, dist = dist) %>%
  t() %>%
  as.data.frame() %>%
  dplyr::rename(functional = 1) %>%
  rownames_to_column(var = "sample") %>%
  mutate(functional = if_else(is.nan(functional), 1, functional))

# Merge all metrics
alpha_div <- richness %>%
  full_join(neutral, by = join_by(sample == sample)) %>%
  full_join(phylogenetic, by = join_by(sample == sample)) %>%
  full_join(functional, by = join_by(sample == sample))

4.3.1 Wild samples

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="0_Wild") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
      ggplot(aes(y = value, x = species, group=species, color=species, fill=species)) +
      geom_boxplot(outlier.shape = NA) +
      geom_jitter(alpha=0.5) +
      scale_color_manual(name="species",
          breaks=c("Podarcis_liolepis","Podarcis_muralis"),
          labels=c("Podarcis liolepis","Podarcis muralis"),
          values=c("#e5bd5b", "#6b7398")) +
      scale_fill_manual(name="species",
          breaks=c("Podarcis_liolepis","Podarcis_muralis"),
          labels=c("Podarcis liolepis","Podarcis muralis"),
          values=c("#e5bd5b50", "#6b739850")) +
      facet_wrap(. ~ metric ) +
      coord_cartesian(xlim = c(1, NA)) +
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      )

4.3.2 Acclimation samples

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="1_Acclimation") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
      ggplot(aes(y = value, x = species, group=species, color=species, fill=species)) +
      geom_boxplot(outlier.shape = NA) +
      geom_jitter(alpha=0.5) +
      scale_color_manual(name="species",
          breaks=c("Podarcis_liolepis","Podarcis_muralis"),
          labels=c("Podarcis liolepis","Podarcis muralis"),
          values=c("#e5bd5b", "#6b7398")) +
      scale_fill_manual(name="species",
          breaks=c("Podarcis_liolepis","Podarcis_muralis"),
          labels=c("Podarcis liolepis","Podarcis muralis"),
          values=c("#e5bd5b50", "#6b739850")) +
      facet_wrap(. ~ metric) +
      coord_cartesian(xlim = c(1, NA)) +
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      )

4.3.3 Antibiotics samples

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="2_Antibiotics") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
      ggplot(aes(y = value, x = species, group=species, color=species, fill=species)) +
      geom_boxplot(outlier.shape = NA) +
      geom_jitter(alpha=0.5) +
      scale_color_manual(name="species",
          breaks=c("Podarcis_liolepis","Podarcis_muralis"),
          labels=c("Podarcis liolepis","Podarcis muralis"),
          values=c("#e5bd5b", "#6b7398")) +
      scale_fill_manual(name="species",
          breaks=c("Podarcis_liolepis","Podarcis_muralis"),
          labels=c("Podarcis liolepis","Podarcis muralis"),
          values=c("#e5bd5b50", "#6b739850")) +
      facet_wrap(. ~ metric) +
      coord_cartesian(xlim = c(1, NA)) +
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      )

4.3.4 Transplant_1 samples

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="3_Transplant1") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
      ggplot(aes(y = value, x = type, group=type, color=type, fill=type)) +
      geom_boxplot(outlier.shape = NA) +
      geom_jitter(alpha=0.5) +
      scale_color_manual(name="type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Control","Hot control", "Treatment"),
          values=c("#76b183","#d57d2c", "#4477AA")) +
      scale_fill_manual(name="type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Control","Hot control", "Treatment"),
          values=c("#76b18350","#d57d2c50", "#4477AA50")) +
      facet_wrap(. ~ metric) +
      coord_cartesian(xlim = c(1, NA)) +
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      )

4.3.5 Transplant_2 samples

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="4_Transplant2") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
      ggplot(aes(y = value, x = type, group=type, color=type, fill=type)) +
      geom_boxplot(outlier.shape = NA) +
      geom_jitter(alpha=0.5) +
      scale_color_manual(name="type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Control","Hot control", "Treatment"),
          values=c("#76b183","#d57d2c", "#4477AA")) +
      scale_fill_manual(name="type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Control","Hot control", "Treatment"),
          values=c("#76b18350","#d57d2c50", "#4477AA50")) +
      facet_wrap(. ~ metric) +
      coord_cartesian(xlim = c(1, NA)) +
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      )

4.3.6 Post-Transplant_1 samples

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="5_Post-FMT1") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
      ggplot(aes(y = value, x = type, group=type, color=type, fill=type)) +
      geom_boxplot(outlier.shape = NA) +
      geom_jitter(alpha=0.5) +
      scale_color_manual(name="type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Control","Hot control", "Treatment"),
          values=c("#76b183","#d57d2c", "#4477AA")) +
      scale_fill_manual(name="type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Control","Hot control", "Treatment"),
          values=c("#76b18350","#d57d2c50", "#4477AA50")) +
      facet_wrap(. ~ metric) +
      coord_cartesian(xlim = c(1, NA)) +
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      )

4.3.7 Post-Transplant_2 samples

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="6_Post-FMT2") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
      ggplot(aes(y = value, x = type, group=type, color=type, fill=type)) +
      geom_boxplot(outlier.shape = NA) +
      geom_jitter(alpha=0.5) +
      scale_color_manual(name="type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Control","Hot control", "Treatment"),
          values=c("#76b183","#d57d2c", "#4477AA")) +
      scale_fill_manual(name="type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Control","Hot control", "Treatment"),
          values=c("#76b18350","#d57d2c50", "#4477AA50")) +
      facet_wrap(. ~ metric) +
      coord_cartesian(xlim = c(1, NA)) +
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      )

4.4 Beta diversity

beta_q0n <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  hillpair(., q = 0)

beta_q1n <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  hillpair(., q = 1)

beta_q1p <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  hillpair(., q = 1, tree = genome_tree)

beta_q1f <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  hillpair(., q = 1, dist = dist)

4.5 Permanovas

4.5.0.1 Load required data

meta <- column_to_rownames(sample_metadata, "Tube_code")

4.5.1 1. Are the wild populations similar?

4.5.1.1 Wild: P.muralis vs P.liolepis

wild <- meta %>%
  filter(time_point == "0_Wild")

# Create a temporary modified version of genome_counts_filt
temp_genome_counts <- transform(genome_counts_filt, row.names = genome_counts_filt$genome)
temp_genome_counts$genome <- NULL

wild.counts <- temp_genome_counts[, which(colnames(temp_genome_counts) %in% rownames(wild))]
identical(sort(colnames(wild.counts)), sort(as.character(rownames(wild))))

wild_nmds <- sample_metadata %>%
  filter(time_point == "0_Wild")

4.5.1.2 Number of samples used

[1] 28

beta_div_richness_wild<-hillpair(data=wild.counts, q=0)
beta_div_neutral_wild<-hillpair(data=wild.counts, q=1)
beta_div_phylo_wild<-hillpair(data=wild.counts, q=1, tree=genome_tree)
beta_div_func_wild<-hillpair(data=wild.counts, q=1, dist=dist)

4.5.1.3 Richness

	Df	SumOfSqs	R2	F	Pr(>F)
species	1	1.631953	0.207198	6.795072	0.001
Residual	26	6.244347	0.792802	NA	NA
Total	27	7.876300	1.000000	NA	NA

4.5.1.4 Neutral

	Df	SumOfSqs	R2	F	Pr(>F)
species	1	2.018797	0.2566342	8.976049	0.001
Residual	26	5.847641	0.7433658	NA	NA
Total	27	7.866438	1.0000000	NA	NA

4.5.1.5 Phylogenetic

	Df	SumOfSqs	R2	F	Pr(>F)
species	1	0.3786052	0.2108638	6.947419	0.001
Residual	26	1.4168908	0.7891362	NA	NA
Total	27	1.7954960	1.0000000	NA	NA

4.5.1.6 Functional

	Df	SumOfSqs	R2	F	Pr(>F)
species	1	0.0787916	0.1594096	4.930642	0.062
Residual	26	0.4154800	0.8405904	NA	NA
Total	27	0.4942716	1.0000000	NA	NA

4.5.2 2. Do the antibiotics work?

4.5.2.1 Acclimation vs antibiotics

treat <- meta  %>% #antibiotics samples giving error when plotting
  filter(time_point == "1_Acclimation" | time_point == "2_Antibiotics")

temp_genome_counts <- transform(genome_counts_filt, row.names = genome_counts_filt$genome)
temp_genome_counts$genome <- NULL

treat.counts <- temp_genome_counts[,which(colnames(temp_genome_counts) %in% rownames(treat))]
identical(sort(colnames(treat.counts)),sort(as.character(rownames(treat))))

treat_nmds <- sample_metadata %>%
  filter(time_point == "1_Acclimation" | time_point == "2_Antibiotics")

4.5.2.2 Number of samples used

[1] 53

beta_div_richness_treat<-hillpair(data=treat.counts, q=0)
beta_div_neutral_treat<-hillpair(data=treat.counts, q=1)
beta_div_phylo_treat<-hillpair(data=treat.counts, q=1, tree=genome_tree)
beta_div_func_treat<-hillpair(data=treat.counts, q=1, dist=dist)

4.5.2.2.1 Richness

	Df	SumOfSqs	R2	F	Pr(>F)
time_point	1	2.033529	0.0957379	6.818959	0.001
species	1	2.327685	0.1095866	7.805342	0.001
individual	26	9.722175	0.4577167	1.253885	0.002
Residual	24	7.157208	0.3369589	NA	NA
Total	52	21.240598	1.0000000	NA	NA

4.5.2.2.2 Neutral

	Df	SumOfSqs	R2	F	Pr(>F)
time_point	1	2.168854	0.1067496	9.686182	0.001
species	1	3.090152	0.1520953	13.800733	0.001
individual	26	9.684304	0.4766554	1.663479	0.001
Residual	24	5.373891	0.2644996	NA	NA
Total	52	20.317201	1.0000000	NA	NA

4.5.2.2.3 Phylogenetic

	Df	SumOfSqs	R2	F	Pr(>F)
time_point	1	2.0671823	0.2208963	21.358337	0.001
species	1	0.8731487	0.0933035	9.021460	0.001
individual	26	4.0949687	0.4375828	1.627293	0.005
Residual	24	2.3228577	0.2482174	NA	NA
Total	52	9.3581574	1.0000000	NA	NA

4.5.2.2.4 Functional

	Df	SumOfSqs	R2	F	Pr(>F)
time_point	1	2.1253769	0.4009905	40.2325710	0.001
species	1	0.0070451	0.0013292	0.1333606	0.777
individual	26	1.9000410	0.3584768	1.3833480	0.182
Residual	24	1.2678545	0.2392035	NA	NA
Total	52	5.3003175	1.0000000	NA	NA

beta_richness_nmds_treat <- beta_div_richness_treat$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(treat_nmds, by = c("sample" = "Tube_code"))


beta_neutral_nmds_treat <- beta_div_neutral_treat$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(treat_nmds, by = c("sample" = "Tube_code"))

beta_phylo_nmds_treat <- beta_div_phylo_treat$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(treat_nmds, by = join_by(sample == Tube_code))

beta_func_nmds_treat <- beta_div_func_treat$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(treat_nmds, by = join_by(sample == Tube_code))

4.5.3 3. Do the FMT work?

4.5.3.1 Comparison between FMT2 vs Post-FMT2

transplant3 <- meta  %>%
  filter(time_point == "4_Transplant2" | time_point == "6_Post-FMT2")

transplant3.counts <- temp_genome_counts[,which(colnames(temp_genome_counts) %in% rownames(transplant3))]
identical(sort(colnames(transplant3.counts)),sort(as.character(rownames(transplant3))))

transplant3_nmds <- sample_metadata %>%
  filter(time_point == "4_Transplant2" | time_point == "6_Post-FMT2")

4.5.3.2 Number of samples used

[1] 53

beta_div_richness_transplant3<-hillpair(data=transplant3.counts, q=0)
beta_div_neutral_transplant3<-hillpair(data=transplant3.counts, q=1)
beta_div_phylo_transplant3<-hillpair(data=transplant3.counts, q=1, tree=genome_tree)
beta_div_func_transplant3<-hillpair(data=transplant3.counts, q=1, dist=dist)

4.5.3.2.1 Richness

	Df	SumOfSqs	R2	F	Pr(>F)
species	1	0.9809525	0.0778442	5.775712	0.001
time_point	1	0.7092107	0.0562799	4.175734	0.001
type	1	1.4479860	0.1149060	8.525540	0.001
individual	25	6.9157236	0.5488022	1.628753	0.001
Residual	15	2.5476146	0.2021678	NA	NA
Total	43	12.6014875	1.0000000	NA	NA

4.5.3.2.2 Neutral

	Df	SumOfSqs	R2	F	Pr(>F)
species	1	1.1101711	0.0895457	7.879642	0.001
time_point	1	0.7363935	0.0593971	5.226687	0.001
type	1	1.9027421	0.1534740	13.505059	0.001
individual	25	6.5351386	0.5271203	1.855374	0.001
Residual	15	2.1133659	0.1704628	NA	NA
Total	43	12.3978112	1.0000000	NA	NA

4.5.3.2.3 Phylogenetic

	Df	SumOfSqs	R2	F	Pr(>F)
species	1	0.1341431	0.1005265	6.921553	0.001
time_point	1	0.1159136	0.0868654	5.980943	0.001
type	1	0.1423573	0.1066822	7.345390	0.001
individual	25	0.6512838	0.4880704	1.344205	0.068
Residual	15	0.2907075	0.2178554	NA	NA
Total	43	1.3344053	1.0000000	NA	NA

4.5.3.2.4 Functional

	Df	SumOfSqs	R2	F	Pr(>F)
species	1	-0.0031096	-0.0030953	-0.1177247	0.826
time_point	1	-0.0045864	-0.0045653	-0.1736359	0.844
type	1	0.0775554	0.0771987	2.9361453	0.145
individual	25	0.5385511	0.5360739	0.8155530	0.617
Residual	15	0.3962105	0.3943880	NA	NA
Total	43	1.0046211	1.0000000	NA	NA

beta_richness_nmds_transplant3 <- beta_div_richness_transplant3$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(transplant3_nmds, by = join_by(sample == Tube_code))

beta_neutral_nmds_transplant3 <- beta_div_neutral_transplant3$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(transplant3_nmds, by = join_by(sample == Tube_code))

beta_phylo_nmds_transplant3 <- beta_div_phylo_transplant3$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(transplant3_nmds, by = join_by(sample == Tube_code))

beta_func_nmds_transplant3 <- beta_div_func_transplant3$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(transplant3_nmds, by = join_by(sample == Tube_code))

p0<-beta_richness_nmds_transplant3 %>%
            group_by(individual) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
        scale_color_manual(values=c("#76b183","#d57d2c", "#4477AA")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y = "NMDS2", x="NMDS1 \n Richness beta diversity") +
                theme_classic() +
                theme(legend.position="top")

p1<-beta_neutral_nmds_transplant3 %>%
            group_by(individual) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
        scale_color_manual(values=c("#76b183","#d57d2c", "#4477AA")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y = "NMDS2", x="NMDS1 \n Neutral beta diversity") +
                theme_classic() +
                theme(legend.position="none")
  
p2<-beta_phylo_nmds_transplant3 %>%
            group_by(individual) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
                scale_color_manual(values=c("#76b183","#d57d2c", "#4477AA")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y= element_blank (), x="NMDS1 \n Phylogenetic beta diversity") +
                theme_classic() +
                theme(legend.position="none")

p3<-beta_func_nmds_transplant3 %>%
            group_by(individual) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
                scale_color_manual(values=c("#76b183","#d57d2c", "#4477AA")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y= element_blank (), x="NMDS1 \n Functional beta diversity") +
                theme_classic()+
                theme(legend.position="none")

4.5.4 4. Are there differences between the control and the treatment group?

4.5.4.1 After 1 week –> Post-FMT1

post1 <- meta %>%
  filter(time_point == "5_Post-FMT1")

post1.counts <- temp_genome_counts[,which(colnames(temp_genome_counts) %in% rownames(post1))]
identical(sort(colnames(post1.counts)),sort(as.character(rownames(post1))))

post1_nmds <- sample_metadata %>%
  filter(time_point == "5_Post-FMT1")

4.5.4.2 Number of samples used

[1] 27

beta_div_richness_post1<-hillpair(data=post1.counts, q=0)
beta_div_neutral_post1<-hillpair(data=post1.counts, q=1)
beta_div_phylo_post1<-hillpair(data=post1.counts, q=1, tree=genome_tree)
beta_div_func_post1<-hillpair(data=post1.counts, q=1, dist=dist)

4.5.4.2.1 Richness

	Df	SumOfSqs	R2	F	Pr(>F)
species	1	0.6488906	0.0757032	2.115638	0.001
type	1	0.5615418	0.0655126	1.830847	0.007
Residual	24	7.3610760	0.8587842	NA	NA
Total	26	8.5715084	1.0000000	NA	NA

4.5.4.2.2 Neutral

	Df	SumOfSqs	R2	F	Pr(>F)
species	1	0.7984743	0.104299	3.065171	0.002
type	1	0.6051778	0.079050	2.323148	0.007
Residual	24	6.2519772	0.816651	NA	NA
Total	26	7.6556293	1.000000	NA	NA

4.5.4.2.3 Phylogenetic

	Df	SumOfSqs	R2	F	Pr(>F)
species	1	0.0700374	0.0635900	1.6594454	0.150
type	1	0.0184254	0.0167292	0.4365646	0.821
Residual	24	1.0129278	0.9196808	NA	NA
Total	26	1.1013906	1.0000000	NA	NA

4.5.4.2.4 Functional

	Df	SumOfSqs	R2	F	Pr(>F)
species	1	-0.0048931	-0.0058792	-0.1628739	0.852
type	1	0.1161466	0.1395538	3.8660898	0.091
Residual	24	0.7210172	0.8663254	NA	NA
Total	26	0.8322707	1.0000000	NA	NA

p0<-beta_richness_nmds_post1 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
        scale_color_manual(values=c("#76b183","#d57d2c", "#4477AA")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y = "NMDS2", x="NMDS1 \n Richness beta diversity") +
                theme_classic() +
                theme(legend.position="top")

p1<-beta_neutral_nmds_post1 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
        scale_color_manual(values=c("#76b183","#d57d2c", "#4477AA")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y = "NMDS2", x="NMDS1 \n Neutral beta diversity") +
                theme_classic() +
                theme(legend.position="none")
  
p2<-beta_phylogenetic_nmds_post1 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
                scale_color_manual(values=c("#76b183","#d57d2c", "#4477AA")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y= element_blank (), x="NMDS1 \n Phylogenetic beta diversity") +
                theme_classic() +
                theme(legend.position="none")

p3<-beta_functional_nmds_post1 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
                scale_color_manual(values=c("#76b183","#d57d2c", "#4477AA")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y= element_blank (), x="NMDS1 \n Functional beta diversity") +
                theme_classic()+
                theme(legend.position="none")

grid.arrange(arrangeGrob(p0,p1,p2,p3))

4.5.4.3 After 2 weeks –>Post-FMT2

post2 <- meta %>%
  filter(time_point == "6_Post-FMT2")

post2.counts <- temp_genome_counts[,which(colnames(temp_genome_counts) %in% rownames(post2))]
identical(sort(colnames(post2.counts)),sort(as.character(rownames(post2))))

post2_nmds <- sample_metadata %>%
  filter(time_point == "6_Post-FMT2")

4.5.4.4 Number of samples used

[1] 28

beta_div_richness_post2<-hillpair(data=post2.counts, q=0)
beta_div_neutral_post2<-hillpair(data=post2.counts, q=1)
beta_div_phylo_post2<-hillpair(data=post2.counts, q=1, tree=genome_tree)
beta_div_func_post2<-hillpair(data=post2.counts, q=1, dist=dist)

4.5.4.4.1 Richness

	Df	SumOfSqs	R2	F	Pr(>F)
species	1	0.8966107	0.1132233	3.515588	0.001
type	1	0.6463814	0.0816245	2.534445	0.002
Residual	25	6.3759661	0.8051521	NA	NA
Total	27	7.9189582	1.0000000	NA	NA

4.5.4.4.2 Neutral

	Df	SumOfSqs	R2	F	Pr(>F)
species	1	0.9398968	0.1224370	4.112304	0.002
type	1	1.0227481	0.1332297	4.474800	0.001
Residual	25	5.7139316	0.7443333	NA	NA
Total	27	7.6765766	1.0000000	NA	NA

4.5.4.4.3 Phylogenetic

	Df	SumOfSqs	R2	F	Pr(>F)
species	1	0.1200820	0.1454923	4.647187	0.002
type	1	0.0592745	0.0718175	2.293931	0.036
Residual	25	0.6459931	0.7826902	NA	NA
Total	27	0.8253496	1.0000000	NA	NA

4.5.4.4.4 Functional

	Df	SumOfSqs	R2	F	Pr(>F)
species	1	0.0116738	0.0168395	0.4229082	0.544
type	1	-0.0085273	-0.0123007	-0.3089208	0.913
Residual	25	0.6900904	0.9954612	NA	NA
Total	27	0.6932368	1.0000000	NA	NA

p0<-beta_richness_nmds_post2 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
        scale_color_manual(values=c("#76b183","#d57d2c", "#4477AA")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y = "NMDS2", x="NMDS1 \n Richness beta diversity") +
                theme_classic() +
                theme(legend.position="top")

p1<-beta_neutral_nmds_post2 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
        scale_color_manual(values=c("#76b183","#d57d2c", "#4477AA")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y = "NMDS2", x="NMDS1 \n Neutral beta diversity") +
                theme_classic() +
                theme(legend.position="none")
  
p2<-beta_phylogenetic_nmds_post2 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
                scale_color_manual(values=c("#76b183","#d57d2c", "#4477AA")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y= element_blank (), x="NMDS1 \n Phylogenetic beta diversity") +
                theme_classic() +
                theme(legend.position="none")

p3<-beta_functional_nmds_post2 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
                scale_color_manual(values=c("#76b183","#d57d2c", "#4477AA")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y= element_blank (), x="NMDS1 \n Functional beta diversity") +
                theme_classic()+
                theme(legend.position="none")

grid.arrange(arrangeGrob(p0,p1,p2,p3))